Spatially Resolved Band Gap and Dielectric Function in Two-Dimensional Materials from Electron Energy Loss Spectroscopy

The electronic properties of two-dimensional (2D) materials depend sensitively on the underlying atomic arrangement down to the monolayer level. Here we present a novel strategy for the determination of the band gap and complex dielectric function in 2D materials achieving a spatial resolution down to a few nanometers. This approach is based on machine learning techniques developed in particle physics and makes possible the automated processing and interpretation of spectral images from electron energy loss spectroscopy (EELS). Individual spectra are classified as a function of the thickness with K-means clustering, and then used to train a deep-learning model of the zero-loss peak background. As a proof of concept we assess the band gap and dielectric function of InSe flakes and polytypic WS2 nanoflowers and correlate these electrical properties with the local thickness. Our flexible approach is generalizable to other nanostructured materials and to higher-dimensional spectroscopies and is made available as a new release of the open-source EELSfitter framework.


Contents
S-2 S1 Pooling and clustering of EELS-SI Let us consider a two-dimensional region of the analysed specimen with dimensions L x × L y where EEL spectra are recorded for n p = n x × n y pixels. Then the information contained within an EELS-SI may be expressed as EELS (E ) , i = 1, . . . , n x , j = 1, . . . , n y , = 1, . . . , n E , (S1.1) where I (i,j) EELS (E ) indicates the recorded total electron energy loss intensity for an energy loss E for a location in the specimen (pixel) labelled by (i, j), and n E is the number of bins that compose each spectrum. The spatial resolution of the EELS-SI in the x and y directions is usually taken to be the same, implying that ∆x = ∆y ≈ L x n x = L y n y . (S1.2) For the specimens analysed in this work we have n p = O(10 4 ) spectra corresponding to a spatial resolution of ∆x ≈ 10 nm. On the one hand, a higher spatial resolution is important to allow the identification and characterisation of localised features within a nanomaterial, such as structural defects, phase boundaries, surfaces or edges. On the other hand, if the resolution ∆x becomes too small the individual spectra become noisy due to limited statistics. Hence, the optimal spatial resolution can be determined from a compromise between these two considerations.
In general it is not known what the optimal spatial resolution should be prior to the STEM-EELS inspection and analysis of a specimen. Therefore, it is convenient to record the spectral image with a high spatial resolution and then, if required, combine subsequently the information on neighbouring pixels by means of a procedure known as pooling or sliding-window averaging.
The idea underlying pooling is that one carries out the following replacement for the entries of the EELS spectral image listed in Eq. (S1.1): EELS (E ) → I By increasing the pooling range d, one combines the local information from a higher number of spectra and thus reduces statistical fluctuations, at the price of some loss on the spatial resolution of the measurement. For instance, d = 3/2 averages the information contained on a 3 × 3 square centered on the pixel (i, j). Given that there is no unique choice for the pooling parameters, one has to verify that the interpretation of the information contained on the spectral images does not depend sensitively on their value. In this work, we consider uniform weights, ω |i −i|,|j −j| = 1, but other options such as Gaussian weights with σ 2 = d 2 as variance are straightforward to implement in EELSfitter. The outcome of this procedure is a a modified spectral map with the same structure as Eq. (S1.1) but now with pooled entries. In this work we typically use d = 3 to tame statistical fluctuations on the recorded spectra.
As indicated by Eq. (S1.1), the total EELS intensity recorded for each pixel of the SI receives contributions from both inelastic scatterings and from the ZLP, where the latter must be subtracted before one can carry out the theoretical interpretation of the low-loss region measurements. Given that the ZLP arises from elastic scatterings with the atoms of the specimen, and that the likelihood of these scatterings increases with the thickness, its contribution will depend sensitively with the local thickness of the specimen. Hence, before one trains the deep-learning model of the ZLP it is necessary to first group individual spectra as a function of their thickness. In this work this is achieved by means of unsupervised machine learning, specifically with the K-means clustering algorithm. Since the actual calculation of the thickness has as prerequisite the ZLP determination, see Eq. (S3.22), it is suitable to use instead the total integrated intensity as a proxy for the local thickness for the clustering procedure. That is, we cluster spectra as a function of which coincides with the sum of the ZLP and inelastic scattering normalisation factors. Eq. (S1.6) is inversely proportional to the local thickness t and therefore represents a suitable replacement in the clustering algorithm. In practice, the integration in Eq. (S1.6) is restricted to the measured region in energy loss.
The starting point of K-means clustering is a dataset composed by n p = n x × n y points, ln N (r) tot , r = 1, . . . , n p , r = i + (n y − 1)j , (S1.7) which we want to group into K separate clusters T k , whose means are given by ln N (k) , k = 1, . . . , K . (S1.8) The cluster means represent the main features of the k-th cluster to which the data points will be assigned in the procedure. Clustering on the logarithm of N (r) tot rather than on its absolute value is found to be more efficient, given that depending on the specimen location the integrated intensity will vary by orders of magnitude.
In K-means clustering, the determination of the cluster means and data point assignments follows from the minimisation of a cost function. This is defined in terms of a distance in specimen S-5 thickness space, given by with d rk being a binary assignment variable, equal to 1 if r belongs to cluster k (d rk = 1 for r ∈ T k ) and zero otherwise, and with the exponent satisfying p > 0. Here we adopt p = 1/2, which reduces the weight of eventual outliers in the calculation of the cluster means, and we verify that results are stable if p = 1 is used instead. Furthermore, since clustering is exclusive, one needs to impose the following sum rule The minimisation of Eq. (S1.9) results in a cluster assignment such that the internal variance is minimised and is carried out by means of a semi-analytical algorithm. This algorithm is iterated until a convergence criterion is achieved, e.g. when the change in the cost function between two iterations is below some threshold. Note that, as opposed to supervised learning, here is it not possible to overfit and eventually one is guaranteed to find the solution that leads to the absolute minimum of the cost function. The end result of the clustering process is that now we can label the information contained in the (pooled) spectral image (for r = i + (n y − 1)j) as follows This cluster assignment makes possible training the ZLP deep-learning model across the complete specimen recorded in the SI accounting for the (potentially large) variations in the local thickness.
The number of clusters K is a free parameter that needs to be fixed taking into consideration how rapidly the local thickness varies within a given specimen. We note that K cannot be too high, else it will not be possible to sample a sufficiently large number of representative spectra from each cluster to construct the prior probability distributions, as required for the Monte Carlo method used in this work. We find that K = 10 for the InSe and K = 5 for the WS 2 specimens are suitable choices. Fig. S1.1 displays the outcome of the K-means clustering procedure applied to the InSe specimen, where each color represents one of the K=10 thickness clusters. It can be S-6 compared with the corresponding thickness map in Fig. 2(d); the qualitative agreement further confirms that the total integrated intensity in each pixel N (i,j) tot represents a suitable proxy for the local specimen thickness.

S-7
S2 A deep-learning model for the zero-loss peak Given that the zero-loss peak background cannot be evaluated from first principles, in this work we deploy supervised machine learning combined with Monte Carlo methods to construct a neural network parameterisation of the ZLP. Within this approach, one can faithfully model the ZLP dependence on both the electron energy loss and on the local specimen thickness. Our approach, first presented in, S1 is here extended to model the thickness dependence and to the simultaneous interpretation of the O(10 4 ) spectra that constitute a typical EELS-SI. One key advantage is the robust estimate of the uncertainties associated to the ZLP modelling and subtraction procedures using the Monte Carlo replica method. S2 The neural network architecture adopted in this work is displayed in Fig. 1(b). It contains two input variables, namely the energy loss E and the logarithm of the integrated intensity ln (N tot ), the latter providing a proxy for the thickness t. Both E and ln (N tot ) are preprocessed and rescaled to lie between 0.1 and 0.9 before given as input to the network. Three hidden layers contain 10, 15, and 5 neurons respectively. The activation state of the output neuron in the last layer, ξ where an exponential function is chosen to facilitate the learning, given that the EELS intensities in the training dataset can vary by orders of magnitude. Sigmoid activation functions are adopted for all layers except for a ReLU in the final layer, to guarantee a positive-definite output of the network and hence of the predicted intensity.
The training of this neural network model for the ZLP is carried out as follows. Assume that the input SI has been classified into K clusters following the procedure of App. S1. The members of each cluster exhibit a similar value of the local thickness. Then one selects at random a representative spectrum from each cluster, S-8 each one characterised by a different total integrated intensity evaluated from Eq. (S1.6), such that (i k , j k ) belongs to the k-th cluster. To ensure that the neural network model accounts only for the energy loss E dependence in the region where the ZLP dominates the recorded spectra, we remove from the training dataset those bins with E ≥ E I,k with E I,k being a model hyperparameter S1 which varies in each thickness cluster. The cost function C ZLP used to train the NN model is then where the total number of energy loss bins that enter the calculation is the sum of bins in each individual spectrum, n E = K k=1 n (k) E . The denominator of Eq. (S2.4) is given by σ k (E k ), which represents the variance within the k-th cluster for a given value of the energy loss E k . This variance is evaluated as the size of the 68% confidence level (CL) interval of the intensities associated to the k-th cluster for a given value of E k .
For such a random choice of representative cluster spectra, Eq. (S2.2), the parameters (weights and thresholds) of the neural network model are obtained from the minimisation of Eq. (S2.4) until a suitable convergence criterion is achieved. Here this training is carried out using stochastic gradient descent (SGD) as implemented in the PyTorch library, S3 specifically by means of the ADAM minimiser. The optimal training length is determined by means of the look-back crossvalidation stopping. In this method, the training data is divided 80%/20% into training and validation subsets, with the best training point given by the absolute minimum of the validation cost function C (val) ZLP evaluated over a sufficiently large number of iterations. In order to estimate and propagate uncertainties associated to the ZLP parametrisation and subtraction procedure, here we adopt a variant of the Monte Carlo replica method S1 benefiting from the high statistics (large number of pixels) provided by an EELS-SI. The starting point is selecting where now the superindices (i m,k , j m,k ) indicate a specific spectrum from the k-th cluster that has been assigned to the m-th replica. Given that these replicas are selected at random, they provide a representation of the underlying probability density in the space of EELS spectra, e.g. those spectra closer to the cluster mean will be represented more frequently in the replica distribution.
By training now a separate model to each of the N rep replicas, one ends up with another Monte Carlo representation, now of the probability density in the space of ZLP parametrisations. This is done by replacing the cost function Eq. (S2.4) by makes possible subtracting the ZLP from the measured EELS spectra following the matching procedure described in S1 and hence isolating the inelastic contribution in each pixel, inel (E) over the MC replica sample estimates the uncertainties associated to the ZLP subtraction procedure. By means of these MC samplings of the probability distributions associated to the ZLP and inelastic components of the recorded spectra, one can evaluate the relevant derived quantities with a faithful error estimate. Note that in our approach error propagation is realised without the need to resort to any approximation, e.g. linear error analysis.
One important benefit of Eq. (S2.6) is that the machine learning model training can be carried S-10 out fully in parallel, rather than sequentially, for each replica. Hence our approach is most efficiently implemented when running on a computer cluster with a large number of CPU (or GPU) nodes, since this configuration maximally exploits the parallelization flexibility of the Monte Carlo replica method.
As mentioned above, the cluster-dependent hyperparameters E I,k ensure that the model is trained only in the energy loss data region where ZLP dominates total intensity. This is illustrated by the scheme of Fig. 3.3 in, S1 which displays a toy simulation of the ZLP and inelastic scattering contributions adding up to the total recorded EELS intensity. The neural network model for the ZLP is then trained on the data corresponding to region I, while region II is obtained entirely from model predictions. To determine the values of E I,k , we evaluate the first derivative of the total recording intensity, dI EELS (E)/dE, for each of the members of the k-th cluster. When this derivative crosses zero, the contribution from I inel will already be dominant. There are then two options. First, one sets E I,k = f × E min,k , where f < 1 and E min,k is the energy where the median of dI EELS /dE crosses zero (first local minimum) for cluster k. Second, one sets E I,k to be the value where at most f % of the models have crossed dI EELS /dE = 0, with f 10%. This choice implies that 90% of the models still exhibit a negative derivative. We have verified that compatible results are obtained with the two choices, indicating that results are reasonably stable with respect to the value of the hyperparameter E I,k .
The second model hyperparameter, denoted by E II,k in Fig. 3.3 in, S1 indicates the region for which the ZLP can be considered as fully negligible. Hence in this region III we impose that I ZLP (E) → 0 by means of the Lagrange multiplier method. This condition fixes the model behaviour in the large energy loss limit, which otherwise would remain unconstrained. Since the ZLP is known to be a steeply-falling function, E II,k should not chosen not too far from E I,k to avoid an excessive interpolation region. In this work we use E II,k = 3 × E I,k , though this choice can be adjusted by the user.
Finally, we mention that the model hyperparameters E I,k and E II,k could eventually be determined by means of an automated hyper-optimisation procedure as proposed in, S4 hence further reducing the need for human-specific input in the whole procedure.

S3 Kramers-Kronig analysis of EEL spectra
Here we provide an overview of the theoretical formalism, based on, S5 adopted to evaluate the single-scattering distribution, local thickness, bandgap energy and type and complex dielectric function from the measured EELS spectra. As indicated by Eq. (1) in the main manuscript, these spectra receive three contributions: the one from inelastic single scatterings off the electrons in the specimen, the one associated to multiple inelastic scatterings, and then the ZLP arising from elastic scatterings and instrumental broadening. Hence a generic EEL spectrum I EELS (E) can be decomposed as where E is the energy loss experienced by the electrons upon traversing the specimen, I ZLP is the ZLP intensity, and I (n) inel indicates the contribution associated to n inelastic scatterings. The ZLP intensity can be further expressed as where R(E) is known as the resolution or instrumental response function whose full width at halfmaximum (FWHM) indicates the resolution of the instrument. The normalisation factor N 0 thus corresponds to the integrated intensity under the zero-loss peak. In the following, we assume that the ZLP contribution to Eq. (S3.1) has already been disentangled from that associated to inelastic scatterings by means of the subtraction procedure described in App. S2.
The single-scattering distribution. If one denotes by t the local thickness of the specimen and by λ the mean free path of the electrons, then assuming that inelastic scatterings are uncorrelated and that t ∼ > λ, one has that the integral over the n-scatterings distribution I with B a normalisation constant. From the combination of Eqns. (S3.1) and (S3.3) it follows that S-12 and hence one finds that the integral over the n-scatterings distribution is such that S3.5) in terms of the normalisation N inel of the full inelastic scattering distribution, the sample thickness t and the mean free path λ. Note also that the ZLP normalisation factor N 0 is then given in terms of the inelastic one as and hence one has the following relations between integrated inelastic scattering intensities In order to evaluate the local thickness of the specimen and the corresponding dielectric function, it is necessary to deconvolute the measured spectra and extract from them the singlescattering distribution (SSD), I SSD (E). The SSD is related to the experimentally measured n = 1 distribution, I inel (E) by the finite resolution of our measurement apparatus: where in the following ⊗ denotes the convolution operation. It can be shown, again treating individual scatterings as uncorrelated, that the experimentally measured n = 2 and n = 3 multiple scattering distributions can be expressed in terms of the SSD as and likewise for n ≥ 4. Combining this information, one observes that the spectrum Eq. (S3.1) can be expressed in terms of the resolution function R, the ZLP normalisation N 0 , and the single-S-13 scattering distribution I SSD as follows where δ(E) is the Dirac delta function. If the ZLP normalisation factor N 0 and resolution function R(E) are known, then one can use Eq. (S3.11) to extract the SSD from the measured spectra by means of a deconvolution procedure.
SSD deconvolution. The structure of Eq. (S3.11) indicates that transforming to Fourier space will lead to an algebraic equation which can then be solved for the SSD. Here we define the Fourier S3.12) whose inverse is given by S3.13) which has the useful property that convolutions such as Eq. (S3.8) are transformed into products, 14) The Fourier transform of Eq. (S3.11) leads to the Taylor expansion of the exponential and hence which can be solved for the Fourier transform of the single scattering distribution By taking the inverse Fourier transform, one obtains the sought-for expression for the single scattering distribution as a function of the electron energy loss with the corresponding inverse transformation being If one approximates the continuous function f (E) by its discretised version f (E 0 + n∆E) and S3.20) and likewise for the inverse transform In practice, the EELS spectra considered are characterised by a fine spacing in E and the discrete approximation for the Fourier transform produces results very close to the exact one.
Thickness calculation. Once the SSD has been determined by means of the deconvolution procedure summarised by Eq. (S3.17), it can be used as input in order to evaluate the local S-15 sample thickness t from the experimentally measured spectra. Kramers-Kronig analysis provides the following relation between the thickness t, the ZLP normalisation N 0 , and the single-scattering distribution, , (S3.22) where we have assumed that the effects of surface scatterings can be neglected. In Eq. (S3.22), a 0 = 0.0529 nm is Bohr's radius, F is a relativistic correction factor, with E 0 being the incident electron energy, (E) is the complex dielectric function, and θ E is the characteristic angle defined by

S-16
an insulator or semi-conducting material using The complex dielectric function. The dielectric function of a material, also known as permittivity, is a measure of how easy or difficult it is to polarise a dielectric material such an insulator upon the application of an external electric field. In the case of oscillating electric fields such as those that constitute electromagnetic radiation, the dielectric response will have both a real and a complex part and will depend on the oscillation frequency ω, S3.28) which can also be expressed in terms of the energy E =hω of the photons that constitute this electromagnetic radiation, In the vacuum, the real and imaginary parts of the dielectric function reduce to Re [ (E)] = 1 and Im [ (E)] = 0. Furthermore, the dielectric function is related to the susceptibility χ by S3.30) where ν is the so-called Coulomb matrix.
The single scattering distribution I SSD (E) is related to the imaginary part of the complex dielectric function (E) by means the following relation S3.31) in terms of the sample thickness t, the ZLP normalisation N 0 , and the microscope operation parameters defined in Sect. S3. We can invert this relation to obtain Since the prefactor in Eq. (S3.32) does not depend on the energy loss E, we see that Im[−1/ (E)] will be proportional to the single scattering distribution I SSD (E) with a denominator that decreases with the energy (since θ E ∝ E) and hence weights more higher energy losses.
Given that the dielectric response function is causal, the real part of the dielectric function can be obtained from the imaginary one by using a Kramers-Kronig relation of the form where P stands for Cauchy's prescription to evaluate the principal part of the integral. A particularly important application of this relation is the E = 0 case, implies that and hence one can express the dielectric function in terms of experimentally accessible quantities, Once the complex dielectric function of a material has been determined, it is possible to evaluate related quantities that also provide information about the opto-electronic properties of a material.
One example of this would be the optical absorption coefficient, given by which represents a measure of how far light of a given wavelength λ = hc/E can penetrate into a material before it is fully extinguished via absorption processes. Furthermore, combining Eqns. (S3.25) and (S3.36) one has that for a semiconductor material, such as those considered in this work, the refractive index is given by the relation which implies a positive, non-zero value of the real part of the complex dielectric function at E = 0.
The complex dielectric function (E) provides direct information on the opto-electronic properties of a material, for example those associated to plasmonic resonances. Specifically, a collective plasmonic excitation should be indicated by the condition that the real part of the dielectric function crosses the x axis, 1 (E) = 0, with a positive slope. These plasmonic excitations typically are also translated by a well-defined peak in the energy loss spectra. Hence, verifying that a plasmonic transition indicated by 1 (E) = 0 corresponds to specific energy-loss features provides a valuable handle to pinpoint the nature of local electronic excitations present in the analysed specimen.
The role of surface scatterings. The previous derivations assume that the specimen is thick enough such that the bulk of the measured energy loss distributions arises from volume inelastic scatterings, while edge-and surface-specific contributions can be neglected. However, for relatively thin samples with thickness t below a few tens of nm, this approximation is not necessarily suitable.
Assuming a locally flat specimen with two surfaces, in this case Eq. (S3.1) must be generalised to with I S (E) representing the contribution from surface-specific inelastic scattering. This surface contribution can be evaluated in terms of the real 1 and imaginary 2 components of the complex dielectric function, where the electron kinetic energy is T = m e v 2 /2.
The main challenge to evaluate the surface component from Eq. (S3.41) is that it depends on the complex dielectric function (E), which in turn is a function of the single scattering distribution obtained from the deconvolution of I inel (E) obtained assuming that I S (E) vanishes. For not too thin specimens, the best approach is then an iterative procedure, whereby one starts by assuming that I S (E) 0, evaluates (E), and uses it to evaluate a first approximation to I S (E) using Eq. (S3.41). This approximation is then subtracted from Eq. (S3.40) and hence provides a better estimate of the bulk contribution I inel (E). One can then iterate the procedure until some convergence criterion is met. Whether or not this procedure converges will depend on the specimen under consideration, and specifically on the features of the EELS spectra at low energy losses, E ∼ < 10 eV. For the specimens considered in this work, it is found that this iterative procedure to determine the surface contributions converges best provided that the local sample thickness satisfies t ∼ > 20 nm.
Validation. The calculations of the local specimen thickness, Eq. (S3.22), and of the complex dielectric function, Eq. (S3.37) presented in this work have been benchmarked with the corresponding implementation available within the HyperSpy framework. S6 Provided one inputs the same inelastic spectra and ZLP parametrisation, agreement between the two calculations is obtained.
This benchmark is illustrated in Fig. S3.1, which compares the EELSfitter-based results with those available from HyperSpy separately for the real and imaginary components of the dielectric function. Both calculations use for the same input ZLP and inelastic spectra, associated to a representative pixel of the WS 2 nanoflower specimen. Residual differences can be attributed to implementation differences e.g. for the discrete Fourier transforms. This validation test further confirms the robustness of the calculations presented in this work.
S-20  Figure S3.1: Comparison of the EELSfitter-based results for the real (left) and imaginary (right panel) components of the dielectric function with those available from HyperSpy for the same input ZLP and inelastic spectrum, associated to a representative pixel of the WS 2 nanoflower specimen.

S-21
S4 Band gap analysis of the EELS low-loss region One important application of ZLP-subtracted EELS spectra is the determination of the bandgap energy and type (direct or indirect) in semiconductor materials. The reason is that the onset of the inelastic scattering intensity provides information on the value of the bandgap energy E bg , while its shape for E ∼ > E bg is determined by the underlying band structure. Different approaches have been put forward to evaluate E bg from subtracted EEL spectra, such as by means of the inflection point of the rising intensity or a linear fit to the maximum positive slope, S7 see also. S8 Following, S1 here we adopt the method of, S9,S10 where the behaviour of I inel (E) in the region close to the onset of the inelastic scatterings is described by  ceases to be valid, as one can also see directly from these plots.

S5 Structural characterisation of the InSe specimen
Here we provide details on the structural characterisation of the n-doped InSe specimens. Each specimen is composed by a InSe nanosheet exhibiting a range of thicknesses. The electronic properties of InSe, such as the band gap value and type, are sensitive to both the layer stacking (β, γ, or ε-phase) as well as to the magnitude and type of doping. S11-S14 In particular, n-doped ε-phase InSe has been reported to exhibit a direct bandgap with value E bg = 1.25 eV. S15 These InSe specimens have been grown by means of the Bridgman-Stockbarger method. Doping with Sn impurities is used to obtain n-type InSe. InSe flakes are obtained from bulk material by the sonication procedure, S16 whereby single InSe crystals are pulverized and added to IPA with a ratio of 2:1 (mg:ml). This combination is then sonicated in a sonic bath for 6 hours while keeping the temperature in the range between 25 • C and 35 • C. The ultra high frequencies lead to gas formation between the layers of the material, building up pressure until adjacent layers split apart.
The flakes in the resulting suspension are then collected and dispersed on top of a TEM grid by pipetting. The structural analysis presented in Fig. S5.2 indicates that the n-doped InSe specimens considered in this work exhibit a crystalline structure characterised by a pure ε-phase. In order to further elucidate the type of the bandgap exhibited by this material, photoluminescence (PL) measurements are carried out. The results, displayed in Fig. S5.3, exhibit a well-defined peak located around 1.26 eV. Hence, we conclude that this material is characterised by a direct bandgap with S-25 energy value E bg ≈ 1.26 eV, consistent with the findings of. S15 Note that PL measurements are characterised by a limited spatial resolution as compared to the STEM-EELS results, and therefore this bandgap value corresponds to an average across the specimen. Hence, PL results are not sensitive to spatially-resolved features in the bandgap map such as those reported in Fig. 3(b) of the main manuscript.
S-26 InSe Photoluminesence Figure S5.3: Photoluminescence spectrum acquired in the n-doped InSe specimen. A well-defined peak, is observed indicating that this material is characterised by a direct bandgap with energy value E bg ≈ 1.26 eV.

S6 Bandgap analysis of 2H/3R WS 2 nanoflowers
Here we apply our new approach to the bandgap analysis of same WS 2 specimen considered in the original study. S1,S17 This specimen consisted on a horizontally-standing WS 2 flake belonging to flower-like nanostructures characterised by a mixed 2H/3R polytypism. While S1 restricted its bandgap analysis to a small subset of individual EELS spectra, here we extend it to the whole specimen and as a byproduct also we provide the local thickness map. The goal is to demonstrate how our updated analysis is consistent with the results presented in. S1 The corresponding results for the dielectric function are presented in App. S7.  the magnitude of the 68% CL interval (corresponding to one standard deviation for a Gaussian distribution) from the Monte Carlo replica sample for each pixel of the SI. One finds that the typical uncertainties δE bg range between 15% and 25%. Finally, Fig. S6.2(c) indicates the lower limit of the 68% CL interval for E bg .
Ref. S1 reported a value of the bandgap of 2H/3R polytypic WS 2 of E bg = (1.6 ± 0.3) eV with a exponent of b = 1.3 +0.3 −0.7 extracted from single EELS spectrum. From the spatially-resolved bandgap maps of Fig. S6.2, one observes how our updated results are in agreement with those from the previous study within uncertainties. Furthermore, this spatially-resolved determination of E bg is in agreement within uncertainties with first-principles calculations based on Density Functional Theory (DFT) of the band structure of 2H/3R polytypic WS 2 . S18 These DFT calculations, which also account for spin-orbit coupling effects, find values of E bg in the range between 1.40 eV and 1.48 eV depending on the settings of the calculation. The DFT predictions are hence consistent with the 68% CL interval for the E bg for a wide region of the specimen, as indicated by Fig. S6.2(c).
Furthermore, inspection of the thickness and bandgap maps, Figs. S6.1(b) and S6.2(a) respectively, reveals an apparent dependence of the value of E bg on the local specimen thickness.
Specifically, the bandgap energies tend to increase in the thinner region of the specimen, with t ≈ 25 nm, and then to decrease as one moves towards the thicker regions with t ≈ 50 nm. While this dependence with the thickness is suggestive of the known property of WS 2 that E bg increases S-29 d)  (b) Same as (a) for the corresponding relative uncertainty δE bg , evaluated as half of the 68% CL interval from the Monte Carlo replica sample. (c) Same as (a) indicating the lower limit of the 68% CL interval for E bg . In the three panels, the maps have been filtered such that only those pixels corresponding to the WS 2 specimen are retained.
when going from bulk to monolayer form, uncertainties remain too large to be able to assign significance to this effect.

S-30
S7 Dielectric function in 2H/3R WS 2 nanoflowers App. S6 characterizes the local thickness and the bandgap energy of the 2H/3R WS 2 nanoflower specimen from S1,S17 across the whole EELS-SI. We now present the corresponding results for the spatially-resolved determination of the real, 1 (E), and imaginary, 2 (E), parts of its complex dielectric function. Fig. S7.1 displays 1 (E) and 2 (E) corresponding to two representative spectra of this WS 2 nanoflower specimen. In this analysis we account for the effects of the surface contributions and the error bands quantify the uncertainties associated to the ZLP subtraction procedure. dielectric function exhibits a crossing, 1 (E) = 0 with a positive slope. These crossings can be interpreted as indicating a phase transition involving a collective electronic excitation, such as a plasmonic resonance. Here we define that a crossing takes place wherever 1 (E) = 0 at the 90% CL as estimated from the Monte Carlo representation. For the selected spectra displayed in Fig. S7.1(a,b), this condition is satisfied for E 22 eV and E 18 eV respectively. These values are consistent with the bulk and surface plasmonic resonances in 2H/3R polytypic WS 2 identified in. S17 for those pixels with ≥ 1 crossings. For the upper region of the specimen (characterised by smaller thicknesses), the first crossing is found to be in the low-loss region, while in the bottom (thicker) region, one has a first crossing at E c 21 eV consistent with the WS 2 bulk plasmon peak. We note that one could also show the values of E c for the subsequent crossings, for those pixels exhibiting more than one crossing.

S-32
S8 Simultaneous determination of (E bg , b) from EELS-SI As discussed in the SI Sect. S4, the ZLP-subtracted EELS spectra can be used to extract the value of the bandgap energy E bg as well as its type (direct or indirect) in semiconductor material by fitting a functional form to the onset region of the subtracted inelastic spectra, where b = 0.5 (1.5) corresponds to a direct (indirect) semiconductor. In the results presented in the main manuscript, in particular in for InSe specimen, motivated by the PL analysis, and b = 1.5 for the WS 2 nanoflower, justified by previous work S1,S17 as well as by independent first-principles DFT calculations. S18 Here we present results of fits were both the bandgap energy E bg and type b are simultaneously determined from the ZLP-subtracted EELS spectra in the case of the InSe specimen.
Already in our original study, S1 it was demonstrated that our method is suitable for the joint extraction of both the bandgap energy E bg and the exponent b simultaneously, but that in general the latter is affected by sizable uncertainties. For instance, for two of the EEL spectra considered there from the WS 2 specimen, we found values b = 1.3 +0.3 −0.7 and b = 1.3 +0.3 −0.4 , consistent with the theoretical expectations for an indirect semiconductor. One option to reduce the model uncertainties on the exponent parameter b would be to increase the pooling degree of the pixels in the SI, which would reduce the fluctuations in the low-loss region at the price of a degradation of the achieved spatial resolution. Here instead we show how we can obtain important information about the values of the exponent b in a fully data-driven manner without compromising the spatial resolution. The strategy is to mask away the pixels where the joint determination of (E bg , b) is too noisy, and hence retain only those pixels where the relative uncertainties on the two parameters is below a precision threshold.  Figure S8.1: Upper panels: the median value (left) and 68% CL relative uncertainty (right panel) for the exponent b in joint fits where the values of (E bg , b) are extracted simultaneously from the low-loss region of the subtracted EELS spectra corresponding to the InSe specimen. Only those pixels where the relative uncertainty in the determination of both b and E bg is below the 50% level are retained. Bottom panels: same, now corresponding to the upper (left) and lower (right panel) ranges of the 68% CL interval for b.
uncertainty in the determination of both b and E bg is below the 50% level are retained, and the rest are masked away. One finds that for around 25% of the pixels it is possible to achieve a reasonable precision in these joint fits of (E bg , b).
Inspection of Figs. S8.1 and S8.2 reveals that the EELS data strongly prefers a value of the exponent around b 0.5, consistent with both the PL measurements and with the expectations of this material being a direct semiconductor, while an alternative scenario with b = 1.5 is strongly disfavored. Furthermore, one can verify that for all other pixels in the SI, including the ones masked away, b = 0.5 is contained within the 68% CL uncertainty range. If we average over all the SI pixels displayed in Fig. S8.1 we find b = 0.50 ± 0.26, confirming that indeed the specimen is a direct semiconductor. In addition, by comparing with Fig. 3 (b,c) of the main manuscript, one finds that the extracted values of E bg are stable irrespectively of whether the exponent b is kept fixed or instead is fitted.
The analysis presented here demonstrates that our method is also suitable to provide spatiallyresolved information about the value of the bandgap type, confirming the value of b obtained from the spatially-averaged PL data. Further work will investigate how to improve the trade-off between spatial resolution and precision in the joint fits of (E bg , b) from the subtracted EELS spectra. InSe − Bandgap Energy (68% CL lower limit) S-35